ここで、偏差siがなるべく小さくするようなa,bを決定するわけです。
偏差siの総和を最小にするというアイディアは すぐに思いつきますが、 si=y*i-yiですから、 負になることも考えられるので、単純に足すというわけにはいきません。 そこで、siを自乗すれば正になりますから、 自乗和を最小にするという結論になるわけです。これが最小自乗法です。
学校でオームの法則についての実験を行い、
結果をレポートにまとめるという課題が出ました。
この実験は、電圧を0〜100[V]まで5[V]刻みに変え、
250[Ω]の抵抗に流れる電流を測定するというものです。
当然、理論的にはI=V/Rです。
| 電圧[V] | 電流[mA] | 理論値[mA] | 誤差[%] |
| 0 | 0.0 | 0 | 0.00 |
| 5 | 19.3 | 20 | -3.61 |
| 10 | 38.1 | 40 | -4.70 |
| 15 | 58.8 | 60 | -1.94 |
| 20 | 76.1 | 80 | -4.82 |
| 25 | 102.2 | 100 | 2.24 |
| 30 | 129.0 | 120 | 7.47 |
| 35 | 144.7 | 140 | 3.37 |
| 40 | 156.9 | 160 | -1.96 |
| 45 | 169.6 | 180 | -5.77 |
| 50 | 202.1 | 200 | 1.03 |
| 55 | 210.6 | 220 | -4.25 |
| 60 | 221.6 | 240 | -7.68 |
| 65 | 255.4 | 260 | -1.76 |
| 70 | 298.7 | 280 | 6.69 |
| 75 | 309.2 | 300 | 3.08 |
| 80 | 306.9 | 320 | -4.10 |
| 85 | 335.0 | 340 | -1.47 |
| 90 | 342.0 | 360 | -5.01 |
| 95 | 363.0 | 380 | -4.46 |
| 100 | 432.2 | 400 | 8.05 |
上記のグラフはご察しの通り、スプレッドシートで生成していますので、
近似直線を引くぐらい簡単な操作でやってくれます。
そのとき、通常用いられている方法は最小自乗法(least square method)
といいます。さらにそれによって、近似直線の傾きとy切片を決定し、
未知のデータを予測することを回帰分析(regression analyze)といいます。
この直線が回帰直線(regression line)です。
最小自乗法は下の図のように、すでに得られているi番目の実データ(yi)と、
回帰直線によって予測される値(y*i)の差(偏差:deviation)を考えます。
点の数がn個なら当然、i=1,2,3,・・・,nです。
ここで、偏差siがなるべく小さくするようなa,bを決定するわけです。
偏差siの総和を最小にするというアイディアは
すぐに思いつきますが、
si=y*i-yiですから、
負になることも考えられるので、単純に足すというわけにはいきません。
そこで、siを自乗すれば正になりますから、
自乗和を最小にするという結論になるわけです。これが最小自乗法です。
参考までですが、これは、
を考えることになります。
最小化するには、この式をa,bで偏微分した式を
連立して解けば求めることができます。
具体的には
となります。コンピュータを使えばあっという間に計算できますよ。
せっかくですから、LPを駆使して最適なa,bを決定することは
できないのでしょうか？
残念ながら、最小自乗法はそのままでは使えません。
＃ 工夫すればできます。それはこちらで。
今回は別なアプローチをとることにしましょう。
上で偏差を自乗したのは、実は"微分しやすい"という理由であって
他の策がないわけじゃないのです。今回は絶対値
をとる、という方針でやってみましょう。
この方法をロバスト回帰といいます。
絶対値を使うと微分不可能点が出てきたりして、解析的には嫌ですが
LPによる定式化ならまったく問題ありません。
定式化すると以下のように書けることはすぐに分かるでしょう。
「a xi + b」は予測値y*iです。
そこから実データyiを引いたものの絶対値が偏差siという訳です。
目的関数は、偏差の絶対値の合計を最小にすることです。
一見するとこのままでよさそうですが、実は大きな問題があります。
制約条件に絶対値がついてしまっています。
絶対値つきの式は線形ではありません。
GLPKにはabsという関数が用意されていますが、
これはパラメータに用いるためのものであって、
変数に用いることはできませんので注意してください。
というわけで、絶対値を外す必要があるわけです。
とは言っても、言われてしまえば非常に簡単です。
結論から言えば
と書き直せばOKです。簡単でしょ？
それでは、なんで上のように書けるかを説明します。
上の式は絶対値の中身が正になった場合のものです。
仮に負になったならば、不等号が「≦」で、siは正ですから
制約を満たします。また、偏差は最小化ですから大きくなりすぎることはありません。
下の式は「-si」から分かるように、中身が負になった場合です。
同様に、正になった場合、不等号が「≧」ですから制約を満たします。
このように、2つの式を用いることで絶対値を外すことができました。
これをGLPKでモデル化すれば、以下のように書けます。 例によって定式化と1対1対応ですから、ほぼそのまんまですね。
#absreg.mod
param n;
param x{1..n};
param y{1..n};
var a;
var b;
var s{1..n};
minimize deviation: sum{i in 1..n}s[i];
s.t. abs1{i in 1..n}: a*x[i] + b - y[i] <= s[i];
s.t. abs2{i in 1..n}: a*x[i] + b - y[i] >= -s[i];
#absreg.dat param n:=21; param: x y := 1 0 0.0 2 5 19.3 3 10 38.1 4 15 58.8 5 20 76.1 6 25 102.2 7 30 129.0 8 35 144.7 9 40 156.9 10 45 169.6 11 50 202.1 12 55 210.6 13 60 221.6 14 65 255.4 15 70 298.7 16 75 309.2 17 80 306.9 18 85 335.0 19 90 342.0 20 95 363.0 21 100 432.2;
次に出力結果の抜粋です。
Problem: absreg
Rows: 43
Columns: 23
Non-zeros: 145
Status: OPTIMAL
Objective: deviation = 166.78 (MINimum)
No. Column name St Activity Lower bound Upper bound Marginal
------ ------------ -- ------------- ------------- ------------- -------------
1 a B 3.932
2 b B -0.18
ロバスト回帰の結果を見ると、傾き3.932, y切片-0.18として
回帰直線を引けば、誤差が最小になるようです。
ちなみに、最小自乗法で求めた場合は傾き3.995, y切片-1.11
ですから、それほど大きな差は見られません。
もっとわかりやすく、グラフにしてみました。
赤線が最小自乗法、青線がロバスト回帰です。
グラフで見ても、両方良さそうな近似になっていますね。
抵抗は250[Ω]でしたので、誤差0ならば傾きは4になるはずです。
この点から見ても、近似の確かさが実感できるでしょうか？
ちなみに、この例の場合、グラフは原点を通りますので、
y切片は0です。そのため、モデル中の変数bは無くした方が
より実際を表していると言えそうです。
#absreg2.mod
param n;
param x{1..n};
param y{1..n};
var a;
var s{1..n};
minimize deviation: sum{i in 1..n}s[i];
s.t. abs1{i in 1..n}: a*x[i] - y[i] <= s[i];
s.t. abs2{i in 1..n}: a*x[i] - y[i] >= -s[i];
これで解いてみると、傾き3.929になりました。
このように、思ったモデルにすぐに書き直せるのが
LPを使う魅力です。解析的に解く場合では、
モデルが変わるごとに作り直さなくてはなりませんので
手間がかかってしまいますね。
無事、実験結果をプロットして、きれいな近似直線を
引くことができました。レポートは無事通ったようです。
めでたし、めでたし。
最後に1つ注意しておくと、回帰分析はどのようなデータにも
適用できてしまいます。ということは、直線の傾向を持たない
結果に対しても無理矢理当てはめが可能だということです。
それで得られた結果には何の意味もありませんので注意してください。
詳しいことは他のサイトや文献を当たれば
豊富な情報が得られると思いますので参考にしてください。